/***************************
Make the interpolation in 2D (sort of..)
checked using IDL program on 7-14-10

*****************************/

using namespace std;

extern vector<double> tab_bol;
//time is the x-axis
//mass is the y-axis
double interpolate_iso_phot( double time_ind, double mass_ind0, double mass_ind1, 
		     long size_x, double ** flux_matrix, long matrix_index) {

  

  //initialize variables
  long fx, cx, fy0, cy0, fy1, cy1; 
  double fracx, fracy0, fracy1;
  double p1, p2, p3, p4;
  double interp_mass1, interp_mass2;
  double pval;
   
  
  //definte some frequently used variables
  fx=(long) floor(time_ind);
  cx=(long) ceil(time_ind);
  fracx=time_ind-fx;
  fy0=(long) floor(mass_ind0);
  cy0=(long) ceil(mass_ind0);
  fy1=(long) floor(mass_ind1);
  cy1=(long) ceil(mass_ind1);
  fracy0=mass_ind0-fy0;
  fracy1=mass_ind1-fy1;
 
//ISO: now i need to keep track of 2 versions of
//each of the mass intervals because the tables are different
  long sfy0=size_x*fy0;
  long scy0=size_x*cy0;
  long sfy1=size_x*fy1;
  long scy1=size_x*cy1;
  
 
  //find the values of the 4 points around our
  //value... specifically p1=upper left, p2=upper right, 
  //p3=lower right, p4=lower left

  p1=flux_matrix[matrix_index][fx+scy0];
  p2=flux_matrix[matrix_index][cx+scy1];
  p3=flux_matrix[matrix_index][cx+sfy1];
  p4=flux_matrix[matrix_index][fx+sfy0];
//	cout<<"--"<<fx<<"\t"<<cx<<"\t"<<scy<<"\t"<<sfy<<endl;
// 	cout<<p1<<"\t"<<p2<<"\t"<<p3<<"\t"<<p4<<endl; 
//	cout<<"fracy fracx  "<<fracy<<"\t"<<fracx<<endl;
  //interpolate along the mass direction
  //This is only a local thing that doesn't know about the convention of 2-1 (RdS)
  interp_mass1=(p1-p4)*fracy0+p4;
  interp_mass2=(p2-p3)*fracy1+p3;
//	cout<<"interp_masses  "<<interp_mass1<<"\t"<<interp_mass2<<endl;  
  //now interpolate along the time direction
  pval=(interp_mass2-interp_mass1)*fracx+interp_mass1;
  
// cout<<pval<<endl;
  //now we return the 10^log(variable)
  //pval=pow(10, pval);
  return pval;
}
